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Abstract 

Stable  queueing  systems  consisting  of  two  groups  of  servers,  laving 
exponential  service  times,  placed  in  tandem  and  separated  by  a finite  buffer, 
are  shewn  to  have  a steady-state  probability  vector  of  matrix-gecmetric 
form.  The  queue  is  stable  as  long  as  the  Poisson  arrival  rate  does  not 
exceed  a critical  value,  which  depends  in  a complicated  manner  on  the  ser- 
vice rates,  the  numbers  of  servers  in  each  group,  the  size  of  the  intermediate 
buffer  and  the  unblocking  rule  followed  when  system  becomes  blocked.  The 
critical  input  rate  is  determined  in  a unified  manner. 

For  stable  queues,  it  is  shown  hew  the  stationary  probability  vector 
and  other  important  features  of  the  queue  may  be  computed.  The  essential 
step  in  the  algorithm  is  the  evaluation  of  the  unique  positive  solution  of  a 
quadratic  matrix  equation. 
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1 . Introduction 

The  queueing  model  consisting  of  two  units  in  series  with  a finite 
intermediate  waiting  roan  has  an  extensive  literature,  dating  back  to  1956 
with  the  work  of  G.  C.  Hunt  [9].  The  study  of  blocking  in  two  or  more  units 
in  series  without  intermediate  waiting  spaces  was  initiated  by  B.  Avi- 
Itzhak  and  M.  Yadin  [2].  Further  contributions  to  this  model  are  due  to 
N.  U.  Prabhu  [17]  and  A.  B.  Clarke  [3]. 

Models  in  which  there  is  a finite  waiting  roan  between  the  two  units  and 
the  service  times  in  the  first  unit  have  a general  distribution  were  dis- 
cussed by  T.  Suzuki  [20],  M.  F.  Neuts  [11,12]  and  K.  Hildebrand  [7], 
basically  using  transform  methods  which  are  not  readily  conputationally 
implemented.  The  thesis  by  I.  Hannibalsson  [5]  utilizes  a buffer  model  to 
represent  a queue  with  delayed  feedback.  The  second  unit  then  represents  a 
holding  stage  for  those  customers  who  will  rejoin  the  queue  in  front  of 
Unit  I.  In  this  paper  and  also  in  that  by  B.  Wong,  W.  Gif fin  and  R.  L. 

Disney  [22],  the  analysis  of  finite  capacity  buffer  models  is  carried  out 
by  fairly  involved  spectral  decompositions  of  the  transition  probability 
matrices.  Related  models,  with  finite  total  numbers  of  customers  were  treated 
in  the  papers  by  K.  L.  Arya  [1]  and  0.  P.  Sharma  [19].  These  papers  also  do 
not  have  an  algorithmic  orientation. 

In  recent  years,  there  has  been  a growing  interest  in  the  development 
of  conputational  methods  to  evaluate  the  stationary  probability  vector  and 
related  quantities  for  tandem  queues  vrith  blocking.  This  interest  came 
primarily  frcm  the  recognition  that  these  models  are  useful  in  the  study  of 
the  behavior  of  subsystems  of  canputers.  In  addition  to  detailed  des- 
criptions of  several  computer-related  applications,  A.  G.  Konheim  and  M. 

Reiser  [10]  propose  an  algorithm  for  the  solution  of  a system  consisting  of 
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two  single-server  units  with  exponential  service  time  distributions.  They 
also  allow  feedback  of  seme  departures  fran  the  second  server  to  the  queue 
in  front  of  the  first  unit . In  [18] , these  same  authors  further  considered 
more  elaborate  forms  of  feedback  and  discussed  additional  applications  in 
computer  modeling.  Iterative  numerical  procedures  of  the  Gauss-Seidel 
type,  such  as  proposed  by  F.  S.  Hillier  and  R.  W.  Boling  [8],  may  also  be 
implemented  for  these  models. 

In  addition,  bounds  on  the  blocking  probability  were  investigated  more 
recently  by  F.  G.  Foster  and  H.  G.  Perros  [4].  A particularly  detailed  study 
of  diffusion  approximations  in  tandan  queues  is  due  to  G.  Newell  [15,16], 

The  recent  paper  by  J.  M.  Harrison  [6]  is  also  relevant  in  this  context. 

It  is  the  purpose  of  this  paper  to  show  that  a large  number  of  buffer 
models  with  exponential  servers  may  be  numerically  solved  in  a unified  way. 
The  key  result  identifies  their  stationary  probability  vector  in  a (modi- 
fied] matrix-geometric  form.  Appropriately  partitioning  that  vector  x as 
(xq,  Xp  ...),  we  shew  that 

_i-r+l  , . , 

Xj  = R , for  i r-1, 

where  r is  the  number  of  servers  in  the  first  unit.  The  matrix  R is  the 
unique  positive  solution  to  a matrix  quadratic  equation.  The  spectral 
radius  of  R is  less  than  one . The  r vectors  Xq , . . . jX^^  j are  also  uniquely 
determined . 

The  approach,  vhich  is  used  here,  is  already  implicit  in  the  thesis  of 
V.  Wallace  [21],  but  the  proofs  are  based  on  further  refinements  and  genera- 
lizations given  in  Neuts  [13,14]. 
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Description  of  the  Model 

The  system  consists  of  units  I and  II  and  a finite  intermediate  buffer. 
Unit  I consists  of  r parallel  exponential  servers,  processing  customers  at 
the  same  rate  a.  In  Unit  II,  c parallel  exponential  servers  process  cus- 
tomers at  the  common  rate  g.  Arrivals  to  Unit  I occur  according  to  a homo- 
geneous Poisson  process  of  rate  A.  (See  Fig.  1.) 

The  servers  in  Unit  II  can  be  active  as  long  as  there  are  customers, 
who  have  completed  a pass  through  Unit  I and  are  requesting  their  service. 
There  are  M-c-1  > 0 places  in  the  buffer,  so  that  at  most  M-l  customers  can 
be  either  waiting  in  the  buffer  or  being  processed  by  one  of  the  servers  in 
Unit  II . If  the  number  of  customers  who  have  completed  a pass  through  Unit  I 
but  have  not  cleared  Unit  II  reaches  M,  one  of  the  servers  in  Unit  I becomes 
blocked . 

Depending  on  the  application,  the  blocicing  of  one  or  more  servers  in 
Unit  I my  affect  the  ability  of  the  unblocked  servers  either  to  accept  a 
customer  for  service  or  to  complete  a service  in  course.  We  stall  assume 
that  when  the  number  of  blocked  servers  in  Unit  I reaches  r* , 1 < r4  <r, 
all  unblocked  servers  in  Unit  I also  cease  service . This  situation  will  be 
referred  to  as  full  blocking. 

- - j*- 

In  a partially  blocked  system,  when  a service  completion  in  Unit  II 
occurs,  one  of  the  blocked  servers  of  Unit  I releases  his  customer  into  the 
buffer.  This  server  may  now  again  initiate  a service. 

Next  we  specify  the  unblocking  rule . In  a fully  blocked  system  there  are 
M + r*  - 1 customers  who  have  completed  a pass  through  Unit  I and  are  re- 
questing service  in  Unit  II.  We  define  an  integer  k*,  0 <_  k'-':  <_  M+r5,t-2. 

When  the  number  of  customers,  who  have  not  been  cleared  by  Unit  II,  drops  to 
k*,  all  interrupted  services  in  the  servers  in  Unit  I resume  and  any  free 
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servers  can  again  initiate  services. 

We  shall  allow  a feedback  loop  of  departures  fran  Unit  II  back  to  the 
queue  in  front  of  Unit  I.  With  probability  O'  =1-0,  0 <_  0 < 1,  a cus- 
tomer who  completes  a service  in  Unit  II  leaves  the  system.  Feedback 
occur  s with  probability  0 . 

In  order  to  concentrate  only  on  parameters  which  have  substantial  sig- 
nificance, we  shall  not  discuss  further  extensions  in  which  customers  may 
leave  the  system  from  Unit  I or  may  enter  feedback  loops  from  the  buffer  to 
Unit  I or  from  Unit  II  to  the  buffer.  The  relevant  matrices  which  govern 
such  cases  can  be  constructed  easily;  the  theorems  and  algorithms  discussed 
below  carry  over  routinely. 

We  also  make  the  standard  independence  assumptions.  All  service  times 
and  interarrival  times  are  mutually  independent  random  variables.  From  a 
numerical  viewpoint,  it  is  routine  to  consider  extensions  such  as  the  case 
where  the  rate  of  the  Poisson  arrival  process  depends  on  the  number  of 
blocked  servers  in  Unit  I,  but  in  carder  not  to  add  to  the  number  of  para- 
meters of  the  model  we  shall  not  pursue  this  topic  further. 

Notational  Convention 


The  material  in  this  paper  involves  a large  number  of  Jacobi  matrices, 
whose  detailed  definitions  require  display.  A matrix  such  as 
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will  be  displayed  as 
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2.  The  Structure  of  the  Markov  Process 

Under  the  assumption  of  exponential  service  times  for  the  servers  in 
Units  I and  II,  the  queueing  model  may  be  described  as  a continuous-parameter 
Markov  chain  on  the  state  space  {(i,j),  i >_  0,  0 <_  j <_N},  where  N is  a 
finite  nonnegative  integer.  The  index  i will  denote  the  number  of  customers 
queued  up  or  in  service  in  Unit  I . Such  customers  will  be  called  I-cus- 
tomers . Upon  completion  of  a pass  through  Unit  I,  a customer  becomes  a 
I I -customer . We  note  that  because  of  the  possibility  of  feedback,  a cus- 
tomer may  be  termed  a I-  or  a II-customer  several  times  in  succession  before 
leaving  the  system.  The  role  played  by  the  index  j is  more  complicated  to 
describe  and  will  be  spelled  out  for  the  specific  cases  discussed  belcw. 

In  all  cases,  however,  the  infinitesimal  generator  P of  the  Markov 
chain  will  have  the  structure  of  a block-tridiagonal  matrix  of  the  form 


(1)  P = A 
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block-partitioned  matrix  correspond  to  the  sets  of  states 
{(i,0) , (i,l) , . . . ,(i,N)},  fari>_0. 

We  shall  now  give  the  detailed  definitions  of  these  blocks  for  various 


models  of  increasing  complexity. 
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Model  A 


Unit  I blocks  as  soon  as  there  are  M II-customers  in  tin  system.  All 
c servers  in  Unit  II  are  busy;  there  are  M-c-1  customers  in  the  waiting 
roan  and  one  server  in  Unit  I has  completed  service  of  a customer  who 
cannot  enter  the  waiting  roam.  Unblocking  occurs  as  soon  as  a departure 
from  Unit  II  occurs.  In  terms  of  our  general  description,  Model  A corres- 
ponds to  r*  = 1,  k*  = M-l. 

In  this  case  N = M.  The  matrices  A^,  A^  and  A^  are  given  by 
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and  for  1 < i < r-1. 


Ai0  = 
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For  0 <_  i <_  r-1,  the  matrices  A^  = ^ , and  the  matrices  A^  are  given  by 
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The  asterisks  correspond  to  the  negative  diagonal  entries,  which  are  such 
that  the  row  sums  of  the  matrix  P are  zero. 

Model  B 

This  model  is  as  the  preceding  one,  except  that  full  blocking  occurs  only 
when  r'--,  1 <_  r*  r,  servers  in  Unit  I are  blocked.  The  index  j now  ranges 
fran  0 to  M + r:':-l  and  denotes  the  number  of  II-customers  in  the  system. 
Unblocking  occurs  again  upon  a subsequent  departure  from  Unit  II,  which  cor- 
responds to  the  case  k*  = M + r5’:-2. 

The  blocks  in  the  partitioned  matrix  P are  now  of  order  M + r'-’: . They 
are  obtained  by  augmenting  the  blocks  in  Model  A in  a systematic  manner. 
Specifically 
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where  the  asterisks  stand  respectively  for  the  entries  -X-cB-(r-l)a  , 
-A-c6-(r-r*+l)a,  -A-c0,  chosen  so  that  the  row  sums  of  P are  zero. 


The  matrices  A.  are  given  by 
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for  1 i £ r-1. 

For  0 < i <_  r-1,  the  matrices  Ai2  = A2  and  the  matrices  A^  are  given  by 

• • • 0 0 • . . 0 0 
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II 

The  asterisks  correspond  to  the  negative  diagonal  entries,  which  are 
such  that  the  rcw  sums  in  the  matrix  P are  zero. 

Model  C 

In  this  model,  we  add  further  complexity  to  Model  B by  assuming  that 
when  full  blocking  occurs , Unit  I does  not  become  unblocked  until  the  number 
of  II-custcmers  drops  to  k* . In  most  cases  of  interest,  we  will  have 
c <_  k*  <_  M+r*-2 , and  in  order  to  limit  the  number  of  variants,  we  shall  assume 
that  this  is  the  case.  The  case  k*  = M+r*-2  corresponds  to  Model  B,  so  we 
only  need  to  discuss  the  cases  where  c <_ k*  <_  M+r*-3. 

We  now  consider  the  indices  j : 

0,  1,  ...,  M-l,  M,  ...,  M+r-'-l,  M+r 2,  M+r*-3,  ...,  F*I. 

The  index  values  with  a bar  correspond  to  the  situations  where  the  Unit  I 
is  blocked,  although  fewer  than  M+r*- 1 II-customers  are  in  the  systan. 

The  blocks  in  the  partitioned  matrix  P are  now  matrices  of  order 
2M+2r*-k*-2 . 
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The  matrices  and  A^,  ^ ^or  this  model  ar'e  obtained  by 

adding  M+r*-k*-2  rows  and  columns  to  the  corresponding  matrices  for  Model 
B.  These  rows  and  columns  are  identically  zero. 

The  matrices  and  A^ , 0 <_  i <_  r-1,  are  obtained  by  adding  M+r*-k*-2 
rows  and  columns  to  the  corresponding  matrices  for  Model  B,  and  changing 
the  row  with  index  M+r*-l.  The  diagonal  elements  to  be  added  are  all  equal 
to  -A-cB.  To  the  right  of  these  diagonal  elements  we  add  an  entry  cB9 ' , 
except  for  the  last  row,  where  the  entry  cBG ’ is  placed  in  the  column 
labeled  k*.  In  the  rew  with  index  M+r*-l  the  entry  cB9 ' should  appear 
immediately  to  the  right,  rather  than  to  the  left  of  the  diagonal  entry. 

All  other  added  elements  are  zero. 

The  matrices  A^  and  A^ » 0 < i < r-1,  are  similarly  modified,  with 

entries  A and  cB9  now  playing  the  role  of  the  quantities  -A-cB  and  cB9 ' . 

We  see  that  in  the  present  model,  the  matrices  are  no  longer  Jacobi 
matrices,  but  remain  highly  structured  sparse  matrices.  The  theoretical 
results  in  this  paper  do  not  depend  on  the  detailed  structure  of  the  blocks 
in  the  partitioned  matrix  P,  but  particularly  when  the  order  of  these 
blocks  becanes  large,  their  sparsity  may  be  exploited  to  economize  on  the 
storage  and  processing  time  requirements  of  the  algorithm. 

3.  Quasi-Birth-and-Death  Process 

Consider  an  irreducible  continuous  parameter  Markov  chain  with  state 
space  {(i,j);  i >0,  0 <_  j _<  N } and  infinitesimal  generator  P of  the  form  (1). 

Let  us  denote  by  x the  vector  of  steady-state  probabilities,  associated  to 
P 5 X ^ = = define  the  conservative  stable  matrix  A by  A = 

Aq  + + A^ . We  assume  that  A is  irreducible  and  denote  by  u its  vector 

of  steady-state  probabilities,  i.e.  tt_  A = 0_,  tt_  e = 1.  Each  canponent  of  n 
is  strictly  positive.  In  the  tandem  queue  models  considered  here,  A will 
obviously  be  irreducible. 
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We  partition  x as  x = (Xq,  x^,  — ),  where  each  vector  x^  has  N+l 

canponents.  We  shall  examine  below  the  existence  of  a solution  of  the  form 

i— . . # 

jo  = x^_-j  R for  i > r-1 , where  R has  a spectral  radius  strictly  less 

tlian  one  (sp(R)  < 1).  For  such  a solution  to  exist,  we  must  have  that 

x„  A„ 


^0  “0,1  + Al, 


0 = °’ 


(2) 


-i  Ai,2  + — i+1  Ai+l,l  + — i+2  Ai+2,0  °’ 

h?-2  Ar-2,2  + ^r-l  Ar-l,l  + 

R1_r+1  (A2  + R + R2  Aq)  = 0, 


Vi 


for  0 <_  i <_  r-3 , 


for  i > r-1. 


We  shall  show  that  in  the  positive  recurrent  case,  a strictly  positive 
probability  vector  x of  the  stated  form  exists,  for  which  the  matrix  R is  a 
nonnegative  irreducible  matrix  of  spectral  radius  sp(R)  less  than  one  and 


such  that  A2  + R A^  + R AQ  = 0. 


We  now  lave  to  make  several  technical  assumptions  that  are  satisfied 
for  the  models  we  consider: 

a.  A-.  is  nonsingular.  By  Wallace  [21],  Theorem  3.1,  a sufficient 
condition  is  that  A^  e < 0_  which  means  that  from  any  state  (i,j),  i ^ r, 
it  is  possible  to  move  in  one  step  to  a state  (i+1,  3 * ) or  (i-1,  j').  By 
Wallace  [21],  Lemma  3.4,  A~^  is  a nonpositive  matrix  with  strictly  negative 


diagonal  elements . 

b.  The  matrix  = -A^  A^\  has  at  least  one  nonzero  element  in  each 
row.  A sufficient  condition  is  that  all  diagonal  entries  of  A 2 are  strictly 
positive,  which  means  that  arrivals  can  occur  when  the  system  is  in  any  state 

(i, j) , i 1 r* 

c.  If  we  define  Cq  = -Aq  A^  and  C = Cq  + C2>  we  assume  that  C is 
irreducible . 

2 

The  equation  A?  + R A^  + R Aq=0,  may  now  be  rewritten  in  the  form 
(3)  R = C2  + R2  Cq. 
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Lemma  1 

The  matrix  C lias  a maximal  eigenvalue  equal  to  one  with  corresponding 
left  and  right  eigenvector,  respectively  proportional  to  tt_  and  e. 

Proof : 

a.  »C:-i  (Aq  + A2)  A'1  = i (A1  - A)  A”1  = w, 
since  tt  A = 0_. 

b.  C A^  e = -(Aq  + A2>  e = (A^  - A)  e = A-^  e_, 
since  A e = 0 . 

Define  the  sequence  (R(n),  n >_  0}  of  matrices  as  follows: 

R(0)  = 0, 

R(n+1)  “ C?  + R(n)^  Cq,  for  n >_  0. 

Theorem  1 

2 

If  2L  Aq  - — — A2  — ’ the  e1uation  R = C2  + R Cq  has  a unique  solution 

R for  which  R > 0,  sp(R)  < 1.  This  solution  is  R(n)  and  sp(R)  = 1. 

— JJ.  30 

2 

1 A0  - * — A2  — ’ the  °luation  R = C2  + R C0  has  a unique  solution 

R for  which  R > 0,  sp(R)  < 1.  This  solution,  '*'im  R(n) , is  the  minimal  sol- 

n-*» 

ution  to  the  equation  ard  is  irreducible. 

Remark 

Before  proving  the  theorem,  we  observe  that  assumptions  b.  and  c.  above 
are  needed  to  prove  that  any  matrix  R satisfying  equation  (3)  is  irreducible. 

In  each  of  the  buffer  models  that  we  consider,  one  can  show  easily  that  in 
fact  R(n)  is  a strictly  positive  matrix. 

Proof 

This  theorem  is  proved  by  repeating  almost  verbatim  the  argument  given 
in  [14],  Theorems  1 and  2,  Lemmas  2,  3 and  4.  V/e  only  indicate  the  main 
steps  of  the  proof  here: 

- For  all  n,  j R(n)  < jr,  hence  sp[R(n)]  <_1,  for  all  n. 
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- For  all  n,  R(n+1)  >_R(n),  therefore,  the  sequence  converges 
monotonically  to  a matrix  R such  that 

R = C2  + R2  CQ  >.  0, 

2L  ^ — 2L»  sp(R)  <_  1. 

2 

- Any  matrix  R >_  0,  with  sp(R)  £ 1,  which  satisfies  R = + R C^, 

is  irreducible;  its  maximum  eigenvalue  n is  a solution  of  the  equation 

n = 0(n),  where  for  0 < z<_l,  0(z)  is  the  maximal  eigenvalue  of  the  irre- 

2 

ducible  matrix  C(z)  = + z Cq. 

- The  equation  z = 0(z),  0 < z <_1,  has  the  unique  solution  n = 1» 
iff  e’(l-)  <_1.  If  0 ’ ( 1— ) > 1,  Hat  equation  has  two  solutions: 

n-^  = n < 1 and  n2  = 1*  Moreover,  0 * ( 1— ) = vC'(l)  u,  where  v and  u are 
respectively  the  left  and  right  Perron-Frobenius  eigenvectors  of  C,  nor- 
malized so  that  v e = 1,  _v  u = 1.  However , _v  = n_  and  u = (m  A^  e)  1 A-^e, 
so  that  O'(l-)  = -2  £ A2  e U A^  e)  Clearly  0 ' C 1— ) <_  1,  if  and  only  if 
-2  21  a2  e >_  7t_  A-^  e,  (as  A^  e <_0)  or  equivalently,  < 

2L  A2  — — ^Ai+A2^  — ~ ^ A~Ao ^ ® = 1 Ao  — * 

- The  matrix  R = R(n)  is  the  minimal  solution  to  the  equation, 

n-xx> 

2 

R = C2  + R Cq,  in  the  set  of  nonnegative  matrices  with  sp(R)  <_  1.  This 
matrix  has  as  its  spectral  radius  the  smallest  positive  root  of  z = 6(z), 
and  is  the  only  such  matrix. 

We  now  discuss  the  characterization  of  the  vectors  Xq,  •••>  ^ 

the  case  where  1 A0  e > A2  e . By  repeating  an  argument  given  in  Neuts  [13], 
one  may  show  that  this  is  the  only  case  for  which  the  Markov  process  P is 
positive  recurrent. 


13 


Let  X"  (Xg > , • • • > x.1 ) > and 


A01  A1I 
A10  A20 


r-3,2  r-2,2 

V2,l  Vl,i  + RA0 


r-1,0 


Lemma  2 


F"  is  an  infinitesimal,  generator. 


Proof 

Since  P is  an  infinitesimal  generator  and  R Ag  > 0,  all  off-diagonal 
elements  of  P:;  are  non-negative . 

To  prove  that  P--  e = 0^,  one  needs  only  consider  the  last  N+l  rows  of 
P-,  since  the  other  rows  are  identical  to  rows  of  P.  However 

Vl,0  - + (Ar-l,l+RV  ^ = “A2  £ + R Ag  e 

= -A2  e + R A[)  e + JQ  Ru  (A^RA^R^ ) e 


= R(I-R)-1  (Aq+A1+A2)  e = 0. 


Lemma  3 


Ar-l,l  + R A0  is  irreducible,  then  P*  is  irreducible. 

Note  that  this  condition  is  not  necessary  but  obviously  holds  in  the 
models  we  consider. 


liroof 

P"  1S  irreducible  if  and  only  if  for  all  (i,j)  and  (i*,j') 

(i,i  - 0,1,  ...,  r-1;  j,  j'  = 0,1,  ...,  N)  there  exists  in  P*  a path  from 
(i,j)  to  (i',j').  As  P is  irreducible,  there  exists  in  P a path  for  all 
such  (i,j)  and  (i',j»). 

a.  The  path  in  P involves  only  states  (i",j")  with  0 < i"  < r-1. 

Then  the  same  path  exists  in  P*. 
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b.  The  path  in  P involves  some  state  (i",j")  with  i"  ^ r. 

Let  (r-l,j^)  and  (r-l^)  be  defined  as  follows:  (r-1,3^)  is  the  last 
state  in  the  path  before  passage  for  the  first  time  through  a state  (i*,j*) 
with  i*  >_  r;  (r-1 , j ^ ^ is  the  first  state  in  the  path  after  passage  far  the 
last  time  through  a state  (i+,j+)  with  i+  ^ r.  By  the  structure  of  P,  those 
states  exist  and  in  fact  i*  = i+  = r. 

If  ^ + R Aq  is  irreducible,  there  exists  a path  from  (r-1, j^)  to 

(r-ljjj)  involving  only  states  (r-l,j")  and,  therefore,  there  exists  in  P* 
a path  fran  (i,j)  to  (i',j'). 

Theorem  2 

Under  the  assumption  of  Lemma  3 and  tt  A^  e > £,  let  R > 0 be  the 

2 

minimal  solution  of  R = + R C^.  Let  x*  = (x*,  x*,  ^ be  a 

solution  of  x*  P*  = then  x*  has  components  all  of  the  same  sign. 
Furthermore  x*  may  be  normalized  by 


r-2 

(4)  z x*  e + x*  , (I-R)  e 
v=0  -v  - -r-1 


1. 


The  vector  x = (x^ , x^ , . . . ) with 


(5)  x^  = x$,  for  0 i i i r-1, 

„i-r+l  , i 

= x^_1  R , for  r-1  <_  i, 

is  the  unique,  strictly  positive  steady-state  probability  vector  of  the 
natrix  P. 

The  proof  is  new  obvious. 

Remark 

Since  R is  irreducible  and  sp(R)  < 1,  (I-R)  ^ exists  and  is  strictly 
positive . 
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Corollary  1 

(6)  R e = e. 

Proof 

R = C2  ♦ R2  C0  = -A2  A"1  -R2  % q1, 

hence 

R Ai  e - -A ^ e - R A^  e, 

and 

-R(Aq+A2>  e = -A2  e -R2  AQ  e, 

(I-R)  A2  e = (I-R)  R Aq  e. 

Since  I-R  is  nonsingular,  Formula  (6)  follows. 

Remark 

In  the  tandem  queue  models,  considered  here,  AQ  e is  the  vector  of 
rates  of  departure  frart  Unit  I when  all  servers  are  busy.  A2  e is  the 
vector  of  rates  of  arrival  to  the  queue  in  front  of  Unit  I.  This 
corollary  shows  that  R plays  a role  similar  to  a traffic  coefficient.  In 
numerical  computations , this  relation  serves  usefully  as  an  accuracy  check 
on  the  evaluation  of  R. 

4.  Explicit  Forms  of  the  Equilibrium  Condition 

For  the  specific  versions  A,  B and  C of  the  buffer  model,  the  equili- 
brium condition  jr AQ  c > £ A2  e , may  be  explicitly  written  in  terns  of  the 
parameters  of  the  model.  Although  the  analytic  expressions  of  these  expli- 
cit forms  are  complicated,  they  all  are  of  the  general  form 
(7)  X < (1-8)  r a V 

where  0 < V < 1,  and  V is  a function  of  all  the  parameters  of  the  model, 
except  for  X and  8. 
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Hie  quantity  (1-0)  ra  is  the  critical  input  rate  of  a system  consisting 
only  of  Unit  I with  a feedback  probability  0 . The  entire  right  hand  side  of 
(7)  may  be  interpreted  as  the  critical  input  rate  X*  to  the  system  under 
consideration.  The  dependence  of  X*  on  the  various  parameters  of  ttie  model 
provides  us  with  a readily  accessible  means  of  comparing  the  effects  of  buffer 
size  and  unblocking  rules.  It  must  be  borne  in  mind,  however,  that  queues  for 
wliich  X is  close  to  or  equal  to  X*  will  exhibit  the  typical,  frequently 
undesirable,  long-range  fluctuations  inherent  in  near-critical  queues. 

Theorem  3 

The  vector  tt_  and  the  equilibrium  condition  £ e > tt.  A?  e,  are  given 
by: 

For  Model  A 


■J  1 

for  1 <_  j <_  c, 
for  c <_  j £ M. 


M-l 

(9)  X < (1-Q)  ra  E it.  = (1-0)  ra 

3=0  3 

For  Model  B (2  <_  r*  <.  r) 


for  1 <_  j <_  c, 


for  c <_  j <_  M, 

for  M+l  j < M+r*-l . 


I 


(11)  A < (1-0)  rot 


= (1-0)  rot 


n-1  r’-'-l  f . \ 

Y n.  + Y (l  - i \ttmx.  , 

j=0  J j=l  l r I 

/•  r*-l  ..  . 

"M+r*-!  r ^M+j-lJ 


For  Model  C (1  <_  r*  < r,  c <_  k*  <_  M+r--2) 

We  shall  give  detailed  formulas  only  for  the  most  useful  case  where 
c 1 <_M-1.  The  equations  for  the  other  cases  are  entirely  similar.  We 

denote  by  ffM+ri-:_2’  the  components  of  ti_  corre spending  to  the  indices 

j = k*+l,  M+r--2 . The  explicit  formulas  far  the  components  of  tt  are 

uninspiringly  complicated,  but  their  numerical  values  may  readily  be  computed 
by  solving  the  linear  equations 


(12) 


ra 

3 16  3-1 

_ ra 

3 eg  3-1* 

- ra 

3 ce  j-l  k-+l  ’ 


for  1 <_  j 1 c, 
for  c <_  3 <_  k*, 
for  k5’!+  l^i  <_  M , 


Vj  = 6 ‘ r)  (§)  Vj-1  - W for  1 < j < r*-2, 

rc-o 


V[+r*-l  = ffM+r*-2 


cc  fra^"'1  r"~1 


where 


. . 1'Y^Y  U fi  - ^ v r»'1  G-~)  htr‘ 

h=l  V c6  * v=l  ' r / v=l  ^ r'h=0  Kc*l 

Finally  ttq  is  obtained  frern  the  normalizing  condition  it  e = 1. 

*M-1  r*-l  / A 1 

(13)  A < (1-0)  ra  Y it  . + £ (l  - 1 |ttm..  ..I  . 

3=n  3 j=1  V.  rj  M^-il 

Proof 

We  shall  only  sketch  Hie  proof  for  Model  C.  The  equations  tt_  A = 0 
nvay  be  written  as 


k*+l  c!  IcBi 


-1 

❖ IT 


O’ 


h+r*-l 
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(14)  — l*oi  Tip  + (5ii ^ = 0, 


ra  n._,  - (ra+jfl)  it.  + (j+l)0n.  = 0, 

3 1 J 1 


for  1 < j < c-1, 


r«  ir . - (ra+cB)  it.  + cB  tr.  +6.  , ...  CB  if.  , =0, 

J-l  3 D+l  3,k-  k-+l 

for  c £ j <_  M-l , 

(r"j)a  Vj-1  - [(r-j~1)a  + ce]  V]  + c(3  Vi+I  = °» 

for  0 <_  j <_  r*-3 , 

(r-r*+2)a  *M+rft_3  - |(r-r*+l)a  + cbJ  "M+r*_2  = (for  r*  >_  2) 


(r-r*+l)a  ~ cB  wM+rM.  = 0, 


^M+r*-!  " nM+r*-2 


These  are  clearly  equivalent  to 


V+l* 


(15)  ra  n.  ^ = min(j,c)3  n ^ , for  1 < j < k:’;, 

ra  TTj_^  = cB  n^  + cB  for  k*+l  £ j <_M, 

(r-j )a  = cB  nM+j  + cB  far  1 < 


for  1 < j £r:':-2, 


(r-r*+l)a  it.,,  ...  „ = cB  nM,  ...  .. 

M+r— 2 M+r--l, 


Vi+r'-’-'-l  = 1TH+r*-2  ***  nk*+l* 


Equating  the  expression  recursively  ccmputed  for  n^+^..;  ^ with  , we 

obtain  the  stated  formula  relating  and  n^. 

Hie  inequality  £ e > £ e,  is  equivalent  to 

c M 

X < Z (ran.  - jB0n.)  + Z (ran.  - c80n.) 


j=c+l 


r*-l  p 

•E  r 

j=i  L 


(r-j)  Vj-i  - ce  0 Vj 


c B 0 (M+r*-k*-2 ) , 


and  by  using  the  equations  (14),  we  obtain  Formula  (13). 
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Remarks 

1.  It  is  preferable  not  to  write  the  geometric  sums  in  (8)  and  (10)  in 
closed  forms,  so  that  we  do  not  have  to  write  separate  expressions  for  the 
case  where  ra  = c8  • 


2.  For  r = c = 1,  and  A chosen,  without  loss  of  generality,  to  be  equal  to 
one,  we  obtain  for  Model  A that  the  queue  will  be  stable  if  and  only  if 


(l-e)a 


This  agrees,  after  elementary  manipulations , with  the  conditions  (2)  for 
a i 6,  and  (3)  for  a = 6,  stated  in  Thm  2 of  Konheim  and  Reiser  ([10],  p. 

334) . A minor  correction  is,  hcwever,  needed  in  the  statement  of  that  theorem. 
Condition  (1),  i.e.  (1-0)  min  («,8)  > 1,  is  claimed  to  be  the  equilibrium 
condition  for  the  system  where  M = °°.  As  it  is  implied  by  one  of  the  other 
conditions,  dependiiig  on  whether  a i 8 or  a = 8,  its  inclusion  in  the 
stability  condition  for  finite  M is  clearly  inappropriate. 


5.  The  System  Considered  at  Service  Completions  in  Unit  I 

Upon  considering  the  numbers  of  I-  and  II-custcmers  inmediately  after 
service  completions  in  Unit  I,  we  obtain  a Markov  cliain  with  the  states 
(i,j),  where  i >_  0 and  j = 1,  ...,  M+r*-l.  In  the  interest  of  notational 
simplicity,  we  shall  preserve  the  ear  Tier  state  space,  but  note  that  since 
service  completions  during  full  blocking  in  Unit  I are  impossible,  the  states 
with  j = 0 and  the  additional  states  corresponding  to  full  blocking  are 
ephemeral.  Our  formulas  will  correctly  assign  "steady-state  probabilities" 
equal  to  zero  to  all  such  states  and  it  will  not  be  necessary  to  adjust 
the  dimensions  of  the  matrices  which  are  involved. 
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Theorem  4 

The  stationary  probability  vector  z = (z^,  z,,  . . . ) of  the  embedded 
Markov  chain  at  service  completions  in  Unit  I is  given  by 


2k  = T^lVl,0’  forO  <k<r-l, 

V- -p+  O 

2k  = T Vi  Ao  = T -Vi  R Ao’  for  k -r-1’ 


where  r is  given  by 
(17) 

The  zero  components  in  the  vectors  z,  are  ignored. 
Proof 


r-1 

l 

i=l 


x. 
— i 


A. 

l 


0^  + 


R (I-R)"1  A, 


-1 

e 


Ihe  formulas  (16)  are  readily  obtained  by  a conditioning  argument  for 
the  elementary  probabilities.  The  quantity  i-1dt  is  clearly  the  elementary 
probability  that  a service  completion  occurs  in  (t,t+dt)  and  the  components 

°f  \+l,0  dt  or  ^+1A0  dt  are  elementary  probabilities  of  transitions 

of  the  type  (k+l,j)  (k,  j ' ) . 

Corollary  2 

The  components  with  j=M+v,  0 <_  v < r*-l,  of  the  vector 

1^1 


(18) 


Z = 


k=0  ^ 


• f,  Ai0  + Vl  R(I"R) 


-1  A 
A0  5 


yield  the  stationary  probabilities  that  upon  completion  of  a service  in 
Unit  I,  v + 1 servers  are  blocked  in  Unit  I. 

It  would  satisfy  higher  standards  of  rigor  to  set  up  explicitly  the 
transition  probability  matrix  of  the  embedded  Markov  chain  and  to  verify 
that  z is  indeed  its  invariant  vector.  In  order  to  avoid  introducing  a 
large  amount  of  extra  notation,  we  shall  only  do  this  for  the  case  r = 1. 

In  the  process  we  shall  also  obtain  a different  formula  for  Z,  which  is  also 


21 


easily  implemented  ard  therefore  provides  us  with  an  accuracy  check  in  numer- 
ical computations. 

Completely  elementary  probability  arguments  yield  that  for  r = 1,  the 
transition  probability  matrix  of  the  embedded  chain  is  given  by 
B, 


(19) 


P = 


0 


B1  B2  B3  BH 


A0  A1  A2  A3  A4 


A0  \ A2  A3 


0 0 


A. 


A1  “2 

*0 


where  A = (-A-,1  A0)n  (-An^  An),  and  B = (-An?-  A„)  A , for  n > 0. 
n 12  10  n 01  2 n — 

The  steady-state  equations  z P = z,  will  be  satisfied  by  the  vectors 

defined  in  Formula  (16)  provided  that 


(20) 


-,k+l 


,-l 


k+1 


•h>+l  « fr 


Xq  R A0  - Xq  R aq  ( AQ1  A2)  A^^.  + l Xq  R Aq  \+1_v  Aq, 

v=l 

for  k > 0 . 


Since  x^  = Xq  R (-A^  A~^),  and  using  the  explicit  form  of  the  matrices 


A^,  this  equation  may  be  equivalently  rewritten  as 


Rk+1 

m 


,-l.k+l  k!1  „v+l 


(21)  Xq  15  - (“a2  A^)  a - Z R (-Ag  Af)  (-A2  A^_  ) ’ J AQ  = 0, 


for  k >_  0. 

In  order  to  see  that  the  matrix  in  square  brackets  in  zero,  write  the 
-1  2 -1 

equation  R = (-A2  A^  ) + R (-AQ  A1  ),  k+1  times.  Multiply  the  v-th  equation 
on  the  left  by  Rv  1 and  on  the  right  by  (-A2  A^)k  V+k  and  sum  up. 

Finally  the  expression  for  t is  obtained  fran  the  normalizing  equation 


z e = 1. 


22 


Corollary  3 

In  the  case  r = 1,  the  vector  Z is  also  given  by 


(22) 


Z = (*  Aq  £) 


-1 


TT  A„ 


T ^0  ^A01+A2 


) (A1+A2)  1 Aq  A - (n  AQ  e)  1 AQ  II  A^J  -1  (A^), 


where  II  is  a matrix  with  M + 1 identical  rows  given  by  it  . 
Proof 


Adding  the  steady-state  equations  for  the  matrix  P,  we  obtain 


(23)  Z 


E 


* 'VV'1  \>i 

v-1 


oj  - Zq  ('a01  A 2)  [(A1+A2)  A0]  Zq  -CA1+A2)  ^Aq 


The  matrix  -iA^+A^)  Aq  is  a stochastic  matrix,  whose  first  column  is  zero. 
All  other  elements  are  strictly  positive.  It  has  the  left  invariant  vector 
tt_  Aq  , whose  first  component  is  zero  and  all  other  components  strictly  posi- 
tive. It  now  follows  readily  frcm  the  theory  of  finite  Markov  chains  that 
the  matrix 

(24)  I + (A.j+A2 ) -1  Aq  + UA0e)-1n  Aq  = CA^+A^1  J7  - (uA^)'1  AQ  U A^J  , 

is  nonsingular.  Adding  (w  AQ  e)_1  Z n AQ  = (tt_  AQ  e)-1  n AQ,  to  both  sides 
of  Equation  (23)  and  replacing  Xq  R (-Aq  A~^)  by  Xq,  we  obtain  the  stated 
formula  after  routine  matrix  manipulations. 

We  note  that  the  formula  assigns  the  correct  value  zero  to  the  first 
component  of  Z.  Verifying  that  the  result  so  obtained  agrees  with  Z_ = 
t Xq  R (I-R)  ■*"  Aq,  (Cor.  2)  provides  a partial  check  on  numerical  computa- 
tions . 

6.  Remarks  on  Numerical  Conputations 

The  solution,  presented  here,  lends  itself  to  a ready  numerical  imple- 
mentation. Efficient  programming,  which  takes  the  high  degree  of  sparsity 
of  the  transition  probability  matrix  into  account,  results  in  substantial 
savings  in  memory  storage  and  execution  times.  This  is  particularly  worth- 
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while  when  the  program  is  to  be  used  to  study  the  design  and  control  aspects 
of  the  model.  In  such  studies,  one  or  more  parameters  of  the  model  need  to 
vary  over  a range  of  values,  which  may  require  a substantial  number  of 
executions  of  the  program.  In  view  of  the  complicated  dependence  of  the  model 
on  each  of  its  parameters,  detailed  numerical  studies  appear  to  be  the  only 
way  of  obtaining  the  hard  qualitative  information  needed  in  problems  of 
design  and  optimization.  Such  a numerical  implementation  is  currently  being 
done.  Hie  qualitative  results  derived  from  it  will  be  discussed  elsewhere. 


The  first  step,  after  ascertaining  that  the  queue  is  stable,  is  to 

compute  the  matrix  R.  This  may  be  done  by  successive  substitutions  in  the 
-12-1 

equation  R = -A^  A-^  - R , starting  with  R = 0 . The  relation 

R Aq  e = A^  e,  proved  in  Lemma  2,  serves  as  an  accuracy  check. 

If  r > 1,  the  vectors  x^,  ...,  are  computed  by  solving  the  system 

of  linear  equations,  discussed  in  Theorem  2.  Since  the  number  of  equations  in 
that  system  may  be  very  large,  it  is  again  desirable  to  take  the  special 
structure  of  its  coefficient  matrix  into  account.  Thi‘  may  be  done  as  follows. 
In  the  system 


(25) 


*0  Aoi  + -1  A10  = - 


x , A . „ + x A,  . + x...  A,.,  n - 0 , 
-\i-l  v-1,2  —v  v,l  — ■ o+l  v+1,0  — 


for  1 < v < r-2 , 

^r-2  Ar-2 ,2  + ^r-1  (Ar-l,l+RV  = 

the  matrices  A 0,  0 £ v <_  r-2,  are  clearly  nonsingular,  so  that,  using  all 
v , z 

but  the  first  equation,  we  can  write  the  vectors  x^,  ...,  xr_2’  as 

x^  = 0 < v < r-2,  where  the  matrices  are  readily  canputed.  The 

first  equation  now  yields 


(26) 


Vi  <co  Aoi  * ci  Aio>  = ^ 


which  together  with  the  normalizing  condition 


, g < a-K>-^  = i; 


uniquely  determines  the  vector  x and  hence  also  the  vectors  xn , . . . , 

r-i  -hj 


* '-2* 

If  there  is  need  to  economize  on  memory  storage,  as  when  r and  the 

order  of  the  matrices  are  large,  we  can  avoid  storing  the  matrices  C , as 

they  may  be  evaluated  recursively.  This  can  be  done  using  only  three  arrays 

of  size  N x N and  one  linear  array  of  length  N.  In  the  latter  the  vector 
r-2  ... 

I C e is  accumulated.  This  does  not  significantly  increase  the  processing 
v=0 

time,  as  the  systems  of  equations  = d,  where  d is  a known  vector, 

are  particularly  easy  to  solve. 

The  simplifications,  discussed  above,  are  particularly  striking  when 
0 = 0,  (no  feedback)  as  the  matrices  A^»  0 < v < r-2,  are  then  scalar 
matrices . 

The  remaining  computations  of  the  vector  x,  of  various  moments  and  of 
the  marginal  queue  length  densities,  as  well  as  the  blocking  probabilities 
are  now  entirely  routine. 
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